import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import datetime as dt
# import statsmodels.api as sm
import scipy.optimize as sco
import numpy.random as npr

# 设置中文字体的显示
plt.rcParams["font.sans-serif"] = ["SimHei"]
plt.rcParams["axes.unicode_minus"] = False
data2 = pd.read_csv('data2.csv', index_col="Date")
data2.index = [dt.datetime.strptime(x, "%Y-%m-%d") for x in data2.index]

# 计算对数收益率
log_returns = np.log(data2.pct_change() + 1)
# 蒙特卡罗采样10000个观测点
number = 10000
stock_num = len(log_returns.columns)

# 用现实数据采用蒙特卡罗方法模拟出更多的数据来测量
# 0-1之间均匀分布的函数
weights = npr.rand(number, stock_num)
# print(weights.shape)
# 每期的投资组合权重，这是资产组合的权重，保证每一次的资产组合成分比重相加都是1
weights /= np.sum(weights, axis=1).reshape(number, 1)

# 求散样本均值
# print(log_returns.mean())
# 计算协方差矩阵，这里含有方差数据
# print(log_returns.cov())
# quit()
# 计算年化收益，这里是用期望值

# 组合的收益
prets = np.dot(weights, log_returns.mean()) * 252
# 计算年化风险，这里是用方差，sqrt是开平方的意思，diag是获取对角线的数据
# 各种组合的方差
pvols = np.diag(
    np.sqrt(np.dot(weights, np.dot(log_returns.cov() * 252, weights.T))))
# 简建立一个函数将收益率、波动率、夏普率封装起来
# print(pvols)


def statistics(weights):
    # 这里之所以要乘以252是因为一年中有252个交易日，我们表格中的数据时单日的收益
    ret = np.dot(weights, log_returns.mean()) * 252
    # 计算年化风险，这里是用方差，sqrt是开平方的意思，diag是获取对角线的数据
    vols = np.sqrt(np.dot(weights, np.dot(log_returns.cov() * 252, weights.T)))
    # print(vols)
    # quit()

    return np.array([ret, vols, ret / vols])


# 建立优化目标函数，注意这里有个负号
def min_sharpe(weights):
    return -statistics(weights)[2]


cons = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1})
bnds = tuple((0, 1) for x in range(stock_num))
opts = sco.minimize(min_sharpe,
                    stock_num * [1 / stock_num],
                    method='SLSQP',
                    bounds=bnds,
                    constraints=cons)
opts['x'].round(3)
statistics(opts['x'].round(3))


def min_volatility(weights):
    return statistics(weights)[1]


opts1 = sco.minimize(min_volatility,
                     stock_num * [1 / stock_num],
                     method='SLSQP',
                     bounds=bnds,
                     constraints=cons)
opts1['x'].round(3)
statistics(opts1['x'].round(3))


# 最大化收益组合
def min_return(weights):
    return -statistics(weights)[0]


opts2 = sco.minimize(min_return,
                     stock_num * [1 / stock_num],
                     method='SLSQP',
                     bounds=bnds,
                     constraints=cons)
opts2['x'].round(3)
statistics(opts2['x'].round(3))

trets = np.linspace(0.04, 0.18, 50)
tvols = []
for tret in trets:
    cons1 = ({
        "type": "eq",
        "fun": lambda x: statistics(x)[0] - tret
    }, {
        "type": "eq",
        "fun": lambda x: np.sum(x) - 1
    })
    res = sco.minimize(min_volatility,
                       stock_num * [1 / stock_num],
                       method='SLSQP',
                       bounds=bnds,
                       constraints=cons1)
    tvols.append(res['fun'])
tvols = np.array(tvols)

plt.figure(figsize=(10, 6))
plt.scatter(pvols, prets, c=prets / pvols, marker='o')
plt.scatter(tvols, trets, c=trets / tvols, marker='x', s=100)
plt.plot(statistics(opts['x'])[1],
         statistics(opts['x'])[0],
         '*',
         markersize=15,
         label='最大夏普率组合')
plt.plot(statistics(opts1['x'])[1],
         statistics(opts1['x'])[0],
         '*',
         markersize=15,
         label='最大波动率组合')
plt.plot(statistics(opts2['x'])[1],
         statistics(opts2['x'])[0],
         '*',
         markersize=15,
         label='最大收益组合')
plt.grid(True)
plt.xlabel("预测波动率")
plt.ylabel("预期收益率")
plt.colorbar(label="夏普率")
plt.legend()
plt.show()
